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I. INTRODUCTION 


The fine splitting is a difference between energies of P 3/2 
and Pi /2 states. For hydrogenic systems it can be obtained 
from the Dirac equation, while for many electron systems 
one needs quantum electrodynamic (QED) theory to consis¬ 
tently describe correlations with relativistic effects. The most 
common many-electron Dirac-like methods M are able to 
achieve two significant digits at most, while experimental pre¬ 
cision is about 6 significant digits fiH]. A much more accu¬ 
rate description of light few-electron systems relies on nonrel- 
ativistic version of QED, called NRQED theory. Relativistic, 
retardation, electron self-interaction, and vacuum polarization 
effects can all be accounted for perturbatively by expansion of 
energy levels in powers of the fine structure constant a, 

E{a) = ma'^ +ma'^ +ma^ 6^^'^ +ma^ 4 -... 

where expansion coefficients may include powers of In a. 
Since these expansion coefficients are expressed in terms of 
the first- and second-order matrix elements of some opera¬ 
tors with the nonrelativistic wave function, the accuracy of the 
numerical calculation strongly depends on the quality of this 
function. Eor example, MCHE calculations O-Sl are accurate 
only to three digits because the wave function is a combina¬ 
tion of Slater determinants and does not satisfy the cusp con¬ 
dition. A much more accurate nonrelativistic wave function 
can be obtained by usii^ an explicitly correlated basis such 
as Hylleraas functions iflOl - U^ . However, three-electron inte¬ 
grals with explicitly correlated functions are much more com¬ 
plicated than two-electron ones. Moreover, the required num¬ 
ber of basis functions has to be much larger in order to achieve 
similar accuracy as for two-electron systems. So the extension 
of QED calculations to a three-electron system is not a sim¬ 
ple task. In our recent works |@, [31 we performed complete 
calculations of higher-order ma® and ma’ In a corrections 
to Li and Be+ 2 P 3/2 — 2Pi/2 fine splitting. Here we aim to 
present in more detail the computational methods. 

The fine structure splitting at the leading order is given 
by the expectation value 
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where g is the exact electron g-factor. The mean value in 
Eq. (|2]l, (...) = ($1 ... 1$) is calculated using the wave func¬ 
tion $ from the stationary Schrodinger equation 

(iJ - P)$ = 0 (4) 


with the nonrelativistic Hamiltonian H in the infinite nuclear 
mass limit 
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The Li and Be+ fine structure in the leading order, including 
finite nuclear mass corrections, has been calculated by using 
the Hylleraas functions in Refs. MM- rThe high accuracy 
is achieved by the use of a relatively large number (about 14 
000) of these functions. All matrix elements are expressed in 
terms of standard and extended Hylleraas int egrals, w hich are 
obtained with the help of recursion relations ImUTI] . 

The situation is different with matrix elements of m a® and 
higher-order operators in the Hylleraas basis, where additional 
classes of complicated integrals appear, for which no efficient 
numerical algorithms are known. Other difficulties arise in 
the evaluation of the second-order matrix element with nearly 
singular operators. The Green function, or equivalently the 
sum over pseudo-states, requires large values of nonlinear pa¬ 
rameters. This causes severe problems with the numerical sta¬ 
bility of recursive algorithms with Hylleraas integrals. We 
overcome this problem by the application of another basis set, 
which consists of the explicitly correlated Gaussian functions. 
We have found |[3[3 that the second-order matrix elements 
can be calculated with high precision when nonlinear param¬ 
eters are globally optimized and a large number of Gaussian 
functions is employed. 
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of spin-dependent operators from the Breit-Pauli Hamilto- 


II. HIGHER-ORDER EINE STRUCTURE 

The ma® correction E^^'^ to the fine structure can be ex¬ 
pressed as the sum of the first- and second-order matrix ele- 
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ments with the nonrelativistic wave function, 
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where the Breit-Pauli Hamiltonian is of the form llH 
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The potentially singular second-order quadratic term with 
h)^'> in Eq. O does not contribute to fine splitting, and thus 
is excluded from further consideration, is an effective 

Hamiltonian of order m a®. Following the derivation in Refs. 
Q and iB, can be represented in the following form 
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and is the vector potential at the position of particle a, 
which is produced by all other particles 


= V 


, , ‘^'^ab 

b^a 


§i3 


ab ' ab 


ab 


7 

Pb 


a (gfa X fabY 
2 m rlf, 


( 20 ) 

In order to further improve theoretical predictions, the 
higher-order m contribution is not neglected but instead 
is approximated by the numerically dominating logarithmic 
part. It is obtained from the analogous correction to the he¬ 
lium fine structure |2^ by dropping the cr® terms be¬ 
cause they do not contribute for states with the total electron 
spin S = 1/2, 


4aog = ^<iog) + 2(i?. 
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The neglected higher-order corrections are the nonlogarithmic 
TO term and the finite nuclear mass corrections to the to a® 
contribution. They will limit the accuracy of our theoretical 
predictions for Li and Be+ fine structure. 
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where £a is the static electric field at the position of particle a 


£a = -WaV = -Za'-E + Ya'^ (19) 
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III. TRANSFORMATION OF MATRIX ELEMENTS 

The expectation value of in Eq. (fT^ is transformed 

initially to a form convenient for numerical calculations with 
^P-states 

6 H, = =1(-ZQi+Q2) 

16 

(24) 

dH 2 = ——D 2 

0 
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6 H 4 = ——{Q 3 EQ 4 ) 

(26) 

SH 3 = — [—Z(Q 5 + Qe) + Qt + Qs] + g 

{-ZP 1 +P 2 ) 

+ g (Qg - Qio - Qs) ~ 

(27) 

SHq = — - (Qii + < 512 ) + P 4 

(28) 

SHy = — (Qi4 — Q 15 ) — Y (*5i 7 + Qis) 

(29) 

+ Y {~Q 4 + Qi 9 ~ Q 20 ) 


where and Pi are defined in Table |I1 Additionally, opera¬ 

tors Qi, Q 2 , and Q 4 are transformed into the sum of the sin¬ 
gular D-term with the Dirac-<5 operator and the regular R-part. 
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Matrix elements with D-terms are calculated with Hylleraas, 
while Gaussian functions are used for R-terms, which ensures 
high numerical precision. 

The second-order contribution is split into parts coming 
from intermediate states with specified angular momentum 
and spin. 


where 
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where l 2 ,ix is a projection into doublet or quartet state X, 
respectively. These contributions are also defined in Table [III 
Most of them can be calculated as they stand. Only the non- 


symmetric {Hp'^/{E — H)' H'^^) matrix element needs nu¬ 
merical regularization due to the high singularity of H^\ 
This is done as follows: is transformed to the regular 

form by the following transformations 
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The overall regularized form of is 

= [H^^\ + {Qa,E-H}, (37) 

where 


a a<b 
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The expectation value of the regularized operator is the same 
as that without regularization. What has changed is the 
second-order matrix element 
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These additional Qt operators together with their expectation 
value are presented in Table HIJ 

The last considered term, the ma^ In a correction from 
Eq. (EB, is represented as 


eIH = ln[(Za)-2] 


(41) 

where Di are defined in Table |T] and in Table |II] The 
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Sir V 

+(< 327 ) (Z (Q29} — (Q30} + 2 (Qai))) (44) 

b<a ^ 

= (-Z (<324) + (Q25} - 2 (Q26} 

4 TT \ 

+(<328) (Z {Q29) — (Qao) + 2 (Qai))) (45) 

are expressed in terms of Qi from Table HIH 

IV. SPIN REDUCTION OF MATRIX ELEMENTS 

The wave function T** of the state in a three-electron 
system is of the form 

= -^A[(t)\ri,r2,r:i) [a(l) P{2) -/3(1) a(2)] a(3)] , 

(46) 
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TABLE I: Expectation values of operators for Li and Be’^ 2^Pj states, {Q) = Kj V with the additional prefactor Kj = {1, —1/2} for 
J = 1/2, 3/2, correspondingly. All digits are significant. 


Operator 

Vl; 

i^Be+ 

Ql = I/a pi X Pa = -2 TT Dl - i?l 

-0.695 207 

-14.464 31 

T>1 = I/a *<70 Pa X S^{ra)Pa 

0.097 730 

2.01013 

Rl = Ea*'7aPSPa X :^PaPa 

0.082 895 

1.834 29 

Q 2 = Ea, 67 ^a' 7 «Pa ^ X Pa = -2 TT D 2 - i?2 

-0.502 754 

-11.065 87 

D 2 = Ea.bT^a^T'^T^ S^{rab)Pa 

0.044 668 

0.980 97 

R 2 = Y.a,b^ai^^PaPa X T^Pa pS 

0.222 098 

4.902 28 

QS = Ea.bT^a O'api X P6 = 2 TT D 3 -f i?3 

0.000 421 

0.73715 

773 = Ea,67^a*'7<*Pb >< ^^ir^b)Pb 

0.017 545 

0.369 31 

R 3 = 'Ea,b^a^^^PbPb X :;^Pbpl 

-0.109 834 

-1.583 29 

Qi = J 2 a,bT^a^°'Pb ^ Xpb = 2 tvD 3 + R 4 

0.281 276 

5.67745 

Ri = Ea,67‘a*'7“PbP*> X T^PbPb 

0.171036 

3.35701 

Qb = T,a,b^a '7- % X 

0.161022 

2.122 84 

Qe Ea,bT^a ^ 0 . -fsT^ {Tab ' Pb) 

0.068 423 

0.858 67 

~ Ea, 67 ^ 0 ,CT^a'7a 7^ 7®^ X p/ 

0.189 027 

2.559 31 

Ea, 67 ^a,C 7 ^a (Tac Pc) 

0.052 774 

0.675 94 

A = Ea,6^a(^a X ^b) ^ 

-0.066 977 

-0.904 13 

A Ea,b7‘a.C7‘a('7a X Ob) ^ 3 ^ ^ 3 ^ 

-0.059 905 

-0.82117 

A = Ea,b7^a('7a X CTb)fpi ^ X Pa 

0.102 287 

1.84114 

Qq = Ea.bT^a * '7a pi ^ Pa X Pb 

-0.126 256 

-3.13189 

<3lO = Ea.bT^a * '7a pi X (f/b ' Pb) Pa 

-0.396 739 

-8.005 14 

Qll = Ea,b7^a,C7^b'7a 7 ^ 73 ^ X p/ 

-0.114 547 

-1.44191 

Qi 2 Ea,b 7 ^a,C 7 ^b'7a “7T73~(fbc Pc) 

0.053 650 

0.682 88 

P4 = Ea,b7.a.C7.a('7a X ab) 

0.059 905 

0.82117 

<3i 4 = Ea.bT^a '7a 7^ X Pa 

-0.041132 

0.033 59 

<3i 5 Ea.b 74 a'^a ^- 7 ^ (rab • Pa) 

-0.144617 

-1.287 70 

— Ea,b74a,C54b'7a 777 7 !^ X p/ 

0.171163 

2.362 41 

QlS Ea.bT^a.CT^b '7a ^,3^ ^,3^ (Ab Pa) 

0.065 529 

0.700 97 

<3i 9 = Ea.b74a 1 <7a Pb 7 ^ Pa X Pfc 

-0.224 280 

-3.050 65 

<320 = Ea.b74a*'7aPb 7 t" X (f/b 'Pa )pb 

ah 

-0.506 006 

-8.526 97 


where A denotes antisymmetrization and 0*(ri,r2,r3) is a 
spatial function with Cartesian index i that comes from any of 
the electron coordinates. The normalization we assume is 




(ft,?’2,^ 3 ) 


(49) 


r2, r3)|T’[ci23<?i*(ri,r2,r3)]> = 1 

i i 

(47) 

where V denotes the sum of all permutations of 1,2, and 
3. The and ^P 3/2 wave functions are constructed us¬ 
ing Clebsch-Gordon coefficients. Expectation values with 
these wave functions can be reduced to spatial expecta¬ 
tion values with algebraic prefactor Kj for J = 1/2, 3/2. 
Namely, the first-order matrix elements with auxiliary nota¬ 
tion {iTi/ 2 , 7 T 3 / 2 } take the form 


xa6-gab|4>) = {l,-l/2}(-2e*^'=) ^ 

a/b ab=12,23.31 

r2, r3)\- Qi^)v\c^23 4''" (50) 


where Ckim coefficients are defined in Table IIVI 


($'|0|$) = {1,1} r2, r3)|QT’[ci23 0*(ri,r2,r3)]) 

(48) 


The spin reduction of the second-order matrix elements is 
more complicated. We shall first introduce the following aux- 
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TABLE II: Second-order contributions to Li and Be^ fine splitting X = — Ki^ 2 ) V, the additional prefactor {A' 1/2 ^^ 3 / 2 } is for 

J = 1/2, 3/2, correspondingly. The numerical uncertainties are about 10”“^ 



Contribution 

| 7 fl/ 2 , 7 T 3/2 I 

VLi 

Ifie-I- 



{ 1 , 0 } 

-0.293 49 

-1.0514(2) 

= 


{0,2/3} 

-0.443 91(3) 

-1.625 3 

X2p = 

{mf + H'-Bm 





X^[Hf]r\^ 

{l,-l/ 2 } 

-0.0217(6) 

-1.435(8) 


5X2p 

{l,-l/ 2 } 

-0.086 80 

-2.169 0 



{1,1/4} 

-0.719 6 ( 6 ) 

-5.803(4) 

II 

A 

($|(LfW + hW) ^ ^hW + 






{1/3, 5/6} 

-0.9012(4) 

-3.475 5(6) 



{3,3/10} 

-0.002 31 

-0.033 2 



{- 1 , 1 / 2 } 

0.006 97 

0.102 5 

2^=00 = 


{0,3/2} 

-0.500 75 

-1.885 6(4) 

= 

($|(hW + ffW) ^ (ifW + 





(<l>|77}^’^77^14-) 

{ 2 , 1 } 

-0.733 27(2) 

-2.625 6(2) 



{ 2 , 1 } 

0.000 08 

0.000 9 



{ 2 ,- 1 } 

0.000 00 

0.000 1 



{0,3} 

-0.000 71 

-0.009 6 

Ti = 

(1>|77i"^(^E.<5ArA|4>) 
(-l’|77i"'Ea[<5ArA].|<I>) 

{l,-l/ 2 } 

-0.028 95 

-0.6471 


<5Ti 

{l,-l/ 2 } 

0.007 99 

0.186 0 

T2 = 

(1>|77i"^(^E,<.5Afa6)|4>> 
(<1-1771,"’(^E6<J5Af.A].l<i>) 

{l,-l/ 2 } 

0.00107 

-0.002 5 


5 Y 2 

{l,-l/ 2 } 

-0.003 66 

0.053 5 


iliary functions, 

= T.Qa'P\c(,%<jyi{n,r2,r,) 


Then, the spin-reduced second-order matrix elements are 


(51) 

(52) 


= 'P[ct23{Q\-Q^2)<l^^iri,r2,r3) 


= r 


C123 (<5i3 -Ql%)<l^^iri,r2,r3) 


(53) 

(54) 

_I 


^CTb-(5b|$) = 


{1:0} /m.. 1 




E-H 




(55) 


(i-l £ ^ d. 2^ E s, ■ d.i*) = 1LL2 


'E-H 


(56) 




_0|$\ = _TTi_ 

(E-H)'^' ' 6 \ {E-H)' 


{1:1} Ar. l^P. 


4'’ 
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(57) 
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TABLE III: Expectation values of additional operators arising from reduction of the second-order matrix elements, (Q) = Kj V with the 
additional prefactor Kj = {1, —1/2} for J = 1/2, 3/2, correspondingly. All digits are significant. 



Operator 

14i 

lBe+ 

Q21 = 

Ea,c ^ 

-0.849 430 

-9.552 24 

Q22 = 

'^a,b^a,c if/ ^ 

-1.432170 

-15.223 86 

Q23 = 

'^a,b^a,c if/ ^ 

0.242 656 

3.250 48 

Q24 = 


-0.400 085 

-4.72157 

Q 25 = 

5I^a,Ma,c<d'^“ if^if/ XPa 

-0.766 998 

-8.959 65 

Q 26 = 

Eo,Ma,c<d'^“ T^i^y XP6 

0.159 671 

2.209 81 

Q 27 = 

E — 

5.638 906 

7.898 02 

Q 28 = 

V — 

^a<b r^i, 

2.096 405 

3.233 41 

Q 29 = 

Ea X Pa 

-0.125 946 

-0.969 13 

Qao = 

Eo,Ma ;;:S^ X Po 

-0.224 641 

-1.659 49 

Q 31 = 

Eo.b^^a 7®^ X Pi, 

ah 

0.038 474 

0.360 85 


TABLE IV: Symmetrization coefficients in matrix elements 


(fc, 1 , m) 

Cklm 

A 

^klm 


FI 

^klm 

^klm 


CkL 

P 

^klm 

(1,2,3) 

2 

1 


0 

0 


2 

0 

(1,3,2) 

-1 

-1 


1 

-1 


-1 

1 

(2,1,3) 

2 

1 


0 

0 


2 

0 

(2,3,1) 

-1 

-1 


-1 

1 


-1 

-1 

(3,1,2) 

-1 

1 


1 

-1 


-1 

1 

(3,2,1) 

-1 

-1 


-1 

1 


-1 

-1 


(‘I’lQ 

[E-ny 4^'"“ 

■Qa\^) = 

{l,-l/2} / 
6 \ 

l^Po 

[E-H)' 

jgVfc ^ife\ 


(58) 


a 

{E-H)' 4^ 

■ Qb\^) 

_ {1,1/4}, 
6 ' 

' f 

\ {E - 

^Po A Jmk^,lm\ 

-H)' 1 

(59) 


(4>|^CTa 

e-h = 

{1/3,5/6} 
6 


*‘Po „• Imk 

-H 


(60) 


a b 


(*iE 




l4p 


ab 


- R 


a<b 


E-H 




c<.d 


{3,3/10} 

6 



ijj 

E-H 



( 61 ) 


a b<c 




(63) 


(^1 E-°H E ■ Q&l^) = 


^E-H ^ 


6 


(64) 
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($1 


a<b 




ab 


E-H 


c<.d 




{ 2 , 1 } 




^mnk 




Imn \ 


(65) 


($|^CTa • Qa 

a 


l*£>o 

E-H 


b<c 


{ 2 ,- 1 } 

6 



E-H 


■ kmn^jlmn 
It ^ 


( 66 ) 




a<b 


i ] QlJ |*\ ^ {Q’ ^{ / ,T,fcj» l^fc iT,»ifc\ 

“■b E-H ^ ' 6 \ E-H^^ / 

c<,d 


(67) 


These formulas, including Kj coefficients, have been oh- The resolvent l/(i7 — iJ) for each angular momentum is rep- 
tained with a computer symbolic program. resented in terms of functions with the appropriate Cartesian 

prefactor, as follows 


V. NUMERICAL CALCULATIONS 


(j^So = ^ijkryy 


(69) 


The spatial function (/) in Eq. (l46l l is represented as a linear 
combination of the Hylleraas or the explicitly correlated 
Gaussians functions 112^ 


3-air^-0,21-2-a3r3-Qi2 7-^2-ai3rj3-a23»'23 


p—airi—a2r2 —asra , 


'23 '31 '23 '1 


^n4 „n5 UQ 


( 68 ) 


In the Hylleraas basis we use six sectors with different val¬ 
ues of nonlinear parameters Wi and a maximum value of 
H = ni + no. + 713 - 1-714 + 715 = 12 ; details are presented in 
Refs. lfl^[l3n . In Gaussian basis we use N = 256, 512,1024, 
and 2048 functions with well-optimized nonlinear parameters 
for each basis function separately. The accuracy achieved for 
nonrelativistic energies is about 10“^^ in Hylleraas and 10“^^ 
in Gaussian bases. 

These nonrelativistic wave functions are used in evaluation 
of matrix elements. Most of the Q and P operators in TablesU 
and Hn] are intractable with present algorithms with Hylleraas 
functions due to difficulties with integrals with inverse powers 
of electron distances, but also due to very lengthy expressions 
in terms of Hylleraas integrals. Thus, we calculate them using 
Gaussian functions; however, with some exceptions. There 
are operators Qi, Q 2 , and Q 4 , the expectation value of which 
is very slowly convergent. Namely, the accuracy achieved is 
as low as 10“^ — 10“^ with as many as 2048 well-optimized 
Gaussian functions. So, to avoid loss of numerical accuracy, 
we represent these operators as the sum of the singular D- 
part and the regular i?-part. The singular D-part, numerically 
dominating, is calculated with Hylleraas functions, while the 
regular R-part, free of singularities, is calculated with a Gaus¬ 
sian basis. This leads to significant improvements in accuracy, 
so the numerical uncertainties do not affect theoretical predic¬ 
tions for the fine structure. Numerical results for all first-order 
matrix elements obtained with the largest basis are presented 
in Table HI and Hill The achieved precision is at least 10“^, 
which is one digit better in comparison to second-order ma¬ 
trix elements described in the following. 

The evaluation of second-order matrix elements is much 
more computationally demanding. First of all, they are ob¬ 
tained only in the Gaussian basis, due to its high flexibility. 


^Po 


= r 


(70) 




2 ^ 


(5*7 

T 





(71) 


c 


+ -l{rirtPrirl) 

Xjk 

_ _/ 2 I I I i I I I i I I 

“ ly Va^bE + + ^Ha^bj 

(5^* ■ 

+ dyi + rir\r[) 

-^Kriri + ryy, + ryA) </> ( 72 ) 


where subscripts a, b, and c refer to any of the electrons in¬ 
cluding the same one. Nonlinear parameters for intermediate 
states are extensively optimized for each second-order sym¬ 
metric matrix element. Moreover, one takes all possible rep¬ 
resentations of angular factors for intermediate states in ap¬ 
propriate proportions to ensure the completeness of the basis. 
Most importantly, the number of Gaussian functions for inter¬ 
mediate states is chosen to be sufficiently high to saturate the 
matrix element. Namely, for a given size N of the external 
wave function, we use 3/2 N elements for all Do- and quar¬ 
tet Fo-states, N elements for quartet Po-states, and 1/2 TV for 
Fo-states. Among all matrix elements, the most demanding 
in terms of optimization was that with intermediate states of 
symmetry ^Po, as the external wave function. Here, the ba¬ 
sis set for the resolvent is divided into two sectors. The first 
sector is built of the known basis functions with the nonlinear 
parameters determined in the minimization of E{2'^P). For 
this purpose we took one of the previously generated basis 
sets of 47 of size equal to N/2. The nonlinear parameters of 
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this basis remain fixed during the optimization in order to en¬ 
sure the accurate representation of the states orthogonal to 
The second sector, of size equal to 3/2N or N for the matrix 
element involving Hb or [HaIv, respectively, consists of basis 
functions that undergo optimization. For the asymmetric ma¬ 
trix elements the basis is combined from two corresponding 
symmetric ones. 

The most computationally demanding matrix element was 
the {[HaIt ^/{E — H)' Hb) term, and it has the slowest nu¬ 
merical convergence in the Gaussian basis. Numerical results 
for matrix elements are summarized in Table HJ The achieved 
precision is about 10“^, one digit less than the first-order ma¬ 
trix elements. In addition, we observe significant cancella- 


TABLE V: Summary of m a® contributions to fine splitting. 


Contribution 

Li 

Be+ 


0.293 49 

1.0514(3) 


-0.295 94(2) 

-1.083 5(1) 

X2p^ 

0.735 0(18) 

11.912(24) 

Xip 

-0.423 5(2) 

-1.340 5(3) 


-0.75113(2) 

-2.828 4(6) 


0.733 34(2) 

2.625 7(2) 


-0.002 13 

-0.028 9(1) 

total second order 

0.28916(19) 

10.308(24) 

SHi 

-0.445 2(16) 

-13.160(3) 

SH2 

0.026 31 

0.5778 

SH4 

0.105 63 

2.405 5(1) 

SH5 

0.150 52 

3.186 6 

SHe 

-0.01160 

-0.130 7 

SHt 

-0.02783 

-0.7578(3) 

total first order 

-0.202 1(16) 

-7.879(3) 

total m a® 

0.0871(24) 

2.429(24) 


tions between S' = 1/2 and S = 3/2 intermediate states. 


and between the first- and second-order terms, see Table |V] 
The final numerical result for the m a® contribution in Table 
rvlis relatively quite small. Regarding the ma^ contribution, 
the second-order term Y is numerically dominant, and con¬ 
tributions from Di terms are an order of magnitude smaller. 
Altogether this correction is only three times smaller than the 
m a® contribution, which is certainly not negligible. 


VI. SUMMARY 

We have performed accurate calculations of the fine struc¬ 
ture in Li and Be+ using the nonrelativistic QED approach 
combined with explicitly correlated basis functions. Relativis¬ 
tic and QED corrections are represented in terms of matrix 
elements of effective operators, which are calculated with a 
highly accurate nonrelativistic wave function. Numerical re¬ 
sults are summarized in Table |Vll We observe an agreement 
with the experimental values. However, our result for Li lies 
below, while for Be+ above experiments of lit] and |@1 re¬ 
spectively. As the sign of all corrections is the same for Li 
and Be'*', this may suggest that one of these experiments un¬ 
derestimated its uncertainty. 

The extension of presented computational approach to other 
systems with more electrons is problematic, due to a lack 
of formulas for the four-electron Hylleraas integrals. There¬ 
fore, achieving similar accuracy for the four electron systems 
would be very challenging. 
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Appendix A: Quantum mechanics of three identical particles 


total spin S' = 1/2 

$ = ■^.4[(/)(ri,f2,f3)[a(l)/3(2)-/3(l)a(2)]Q;(3)] 

= ^[a(l) /3(2) a(3) t/i + /3(1) a(2) a(3) '02 + 
v 6 

a(l)a(2)/3(3)03] (A3) 

where A denotes antisymmetrization, and 


01 = 

= 0(01, 

02, 

,03) - 

- 0(02, 

,03, 

,0i) 

+ 0(02, 

,01,03) 

- 


0(03, 

02, 

,0i) 






(A4) 

02 = 

= 0(03, 

01, 

,02) - 

- 0(02, 

,01, 

,03) 

- 0(01, 

,02,03) 

+ 


0(01, 

03, 

,02) 






(A5) 

03 = 

= 0(02, 

03, 

,0i) - 

- 0(03, 

,01, 

,02) 

+ 0(03, 

,02,0l) 

- 


0(01, 

03, 

,02) 






(A6) 


Consider a wave function of three identical particles 
■'Vj Lt)- Eigenstates of the nonrelativistic Hamiltonian 
can be classified by representation of the permutation group 
S 3 . Two of them, id and sgn, are one dimensional, and the 
third is two dimensional. The wave functions corresponding 
to one-dimensional representations are 


0s(Ti,'r2,r3) = ^ [0(^1,02,03) + 0(^2,03, n)-I- (Al) 
V6 

0(03,01,02) + 0(02,01,03) + 0(03, 02, 0l) + 0(01,03,02)] 


and 


0a( 01,02,03) = [0(01,02,03) + 0(02,03, 0l) + (A2) 

0(03,01,02) - 0(02,01,03) - 0(03, 02,0l) - 0(01,03,02)] 


In order to construct the wave functions corresponding to 
the two-dimensional representation, let us consider the spin- 
dependent wave function for a three-electron system for the 


functions form a two-dimensional representation of S 3 , 

Ei 0i = 0. 

Let us denote the standard matrix element 

(0'10)s = (0'X’'i, ^ 2 , r 3 )|T’[ci 23 0*(ri,r2,r3)]) (A7) 

where V denotes the sum of all permutations of 1,2, and 3. 
Then 

($']$) = (0'|0)s (A8) 

and the scalar products between ipi is 

(0il'0i) = (0'l0)s (-1 + 3Sij) (A9) 

The two orthogonal and normalized functions can be chosen 
as 1 /)/ = and 0// = (02 - 03 )/v^- 

The first-order matrix elements of the spin-independent op¬ 
erator Q are 

($|Q|$) = i(0-|Q|00 = (0|Q|0)s, (AlO) 
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and the second-order matrix elements with Qi and Qi are 

= {<^\Qi _ jj Q 2 |</>)s (All) 

In the numerical evaluation of second-order matrix elements 
with doublet S' = 1/2 intermediate states, the resolvent 


1/{E — H) is represented on the basis of functions of proper 
S 3 symmetry, namely ipj and xpn 

(A12) 

= {4>^\E-H\4})s = ENM-nM 

Hence, the second-order matrix element using Eq. (A7) be¬ 
comes 




{<j>\Q,\4>'^)s{EN-'H)l^ (</'|g2|</)s 

(^[ci23 /'(r-i, r 2 , r 3 )|gi {E Af - 'H)/'/ {(/)’’\Q2\'P[ci23 (/(ri, ra, r 3 )) 


and the last form is used in the numerical calculations. 



